

import pandas as pd
import numpy as np
from sklearn.ensemble import IsolationForest  # 异常检测
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# 1. 模拟原始数据集（基于实际业务指标）
data = {
    "timestamp": pd.date_range("2025-07-10 00:00:00", periods=3600, freq="s"),
    "vehicle_id": np.random.choice(["V001", "V002", "V003"], 3600),
    "speed": np.concatenate([np.random.normal(60, 5, 3000), np.random.normal(100, 10, 600)]),  # 混入超速数据
    "acceleration": np.random.uniform(-3.5, 3.5, 3600),  # 急加速/急刹车
    "brake_status": np.random.choice([0, 1], 3600, p=[0.85, 0.15]),  # 制动信号
    "engine_rpm": np.random.randint(800, 3000, 3600),
    "gps_lat": np.cumsum(np.random.normal(0, 0.001, 3600)),  # 模拟行驶轨迹
    "gps_lon": np.cumsum(np.random.normal(0, 0.001, 3600))
}
df = pd.DataFrame(data)

# 2. 数据预处理（Pandas核心操作）
# 异常值清洗：速度>120km/h或加速度绝对值>3.5 m/s²视为异常
df = df[(df["speed"] <= 120) & (df["acceleration"].abs() <= 3.5)]
# 时间序列特征：提取小时与分钟
df["hour"] = df["timestamp"].dt.hour
# 标记夜间驾驶（疲劳驾驶风险）
df["is_night"] = df["hour"].apply(lambda x: 1 if x in [22, 23, 0, 1, 2, 3, 4, 5] else 0)

# 3. 驾驶行为特征工程（Numpy数值计算）
# 急加速：加速度>2.5 m/s²
df["hard_accel"] = np.where(df["acceleration"] > 2.5, 1, 0)
# 急刹车：加速度<-2.5 m/s²
df["hard_brake"] = np.where(df["acceleration"] < -2.5, 1, 0)
# 超速：速度>100km/h
df["overspeed"] = np.where(df["speed"] > 100, 1, 0)

# 4. 聚合驾驶员风险评分（按车辆分组）
driver_risk = df.groupby("vehicle_id").agg(
    total_events=("hard_accel", "count"),  # 总事件数
    risk_ratio=("hard_accel", lambda x: (x.sum() + df.loc[x.index, "hard_brake"].sum() + df.loc[x.index, "overspeed"].sum()) / len(x))
).reset_index()

# 5. UBI风险评级模型（机器学习集成）
# 特征标准化
scaler = StandardScaler()
features = scaler.fit_transform(driver_risk[["risk_ratio", "total_events"]])
# 异常驾驶检测（高风险驾驶员）
model = IsolationForest(contamination=0.1)  # 假设10%高风险
driver_risk["risk_level"] = model.fit_predict(features)
driver_risk["risk_level"] = driver_risk["risk_level"].map({-1: "高风险", 1: "低风险"})

# 6. 可视化与业务输出
plt.figure(figsize=(10, 6))
plt.scatter(driver_risk["total_events"], driver_risk["risk_ratio"], c=[0 if r=="低风险" else 1 for r in driver_risk["risk_level"]], cmap="RdYlGn")
plt.xlabel("总驾驶事件数")
plt.ylabel("高风险事件比例")
plt.title("驾驶员UBI风险评级分布")
plt.colorbar(label="风险等级")
plt.savefig("driver_risk_plot.png")  # 输出至保险报告